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The hydrodynamic equation of a spinor Bose-Einstein condensate (BEG) gives a simple description 
of spin dynamics in the condensate. We introduce the hydrodynamic equation of a ferromagnetic 
BEG with dissipation originating from the energy dissipation of the condensate. The dissipative 
hydrodynamic equation has the same form as an extended Landau-Lifshitz-Gilbert (LLG) equation, 
which describes the magnetization dynamics of conducting ferromagnets in which locahzed magne- 
tization interacts with spin-polarized currents. Employing the dissipative hydrodynamic equation, 
we demonstrate the magnetic domain pattern dynamics of a ferromagnetic BEG in the presence and 
absence of a current of particles, and discuss the effects of the current on domain pattern forma- 
tion. We also discuss the characteristic lengths of domain patterns that have domain walls with and 
without finite magnetization. 

PACS numbers: 03.75.Kk, 03.75.Mn 03.75.Lm 



I. INTRODUCTION 



A particular feature of superfluids and superconductors with spin degrees of freedom, such as superfluid Helium 
three, p-wave superconductors, and spinor Bose-Einstein condensates (BECs) of ultra-cold atoms, is that they support 
the non-dissipative flow of spins, or spin supercurrent [1, 2]. In such systems, the condensed state is described with 
a multi-component order parameter; the supercurrent of the particles in each spin state is proportional to the phase 
gradient of the corresponding component of the order parameter; hence, the gradient of the relative phase of the order 
parameters in different spin states yields the supercurrent of spins. In particular, when the system is spontaneously 
magnetized, the spin supercurrent is expected to give a nontrivial effect on the magnetization dynamics. This is the 
case for a ferromagnetic BEC. In recent experiments, in situ techniques for the imaging of magnetization profiles 
enables us to investigate the real-time dynamics of magnetizations, such as spin texture formation and the nucleation 
of spin vortices, in ferromagnetic BECs [3-5]. 

For the investigation of the magnetization dynamics, a hydrodynamic equation has an advantage. It provides the 
simple description of magnetization dynamics in a ferromagnetic BEC to investigate instabilities [6, 7] and config- 
urations of skyrmions and spin textures [8, 9]. The hydrodynamic equation, in the absence of energy dissipation, 
takes the same form as the Landau-Lifshitz-Gilbert (LLG) equation without damping if the partial time derivative 
is replaced by the material derivative Dt = dt + ^^mass • V, or equivalently, if the adiabatic spin-transfer torque is 
added [10, 11]. Here, Vmass is the superfluid velocity, which is related to the magnetization direction / as 

V X 7;„ass = ^/ • (V/ X V/), (1) 

where F and M are the spin and mass of an atom, respectively. This identity comes from the continuous spin-gauge 
symmetry of the ferromagnetic BEC [12, 13], and is known as the Mermin-Ho relation [14]. Since the spin is transferred 
along the velocity field I'mass [see Eq. (24)], the appearance of I'mass in the hydrodynamic equation is the consequence 
of the spin supercurrent. 

In this paper, we investigate the effect of Vmass on the magnetization dynamics using the hydrodynamic description. 
We show that in the presence of energy dissipation, the Gross-Pitaevskii (GP) equation for a ferromagnetic BEC is 
reduced to a dissipative hydrodynamic equation, which is equivalent to the extended LLG equation written as 

% = lf^ -Beff - T'f X If - (^^a,, . V)/, (2) 

where B^s is an effective magnetic field and F' is a damping parameter. The standard LLG equation, which is 
widely used to describe the magnetization dynamics in ferromagnets [15], consists of a spin torque due to the effective 
magnetic field and a damping term, and it corresponds to Eq. (2) without the third term on the right hand side. 
The extended LLG equation was introduced to describe the magnetization dynamics affected by spin currents [16], 
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which includes additional torque terms, the so-called adiabatic [10, 11] and nonadiabatic [17] spin-transfer torques. 
The third term on the right-hand side of Eq. (2) corresponds to the adiabatic spin torque term. Thus, Eq. (2) is the 
extended LLG equation without the nonadiabatic spin-transfer torque. In conducting ferromagnets, currents can be 
controlled by the external field as well as generated by magnetic texture dynamics [18]. In the case of a ferromagnetic 
BEG, the spin-transfer torques are related to the superfluid velocity fmass) which is induced by spin textures [see 
Eq. (1)]. The damping parameter T' in Eq. (2) corresponds to the so-called Gilbert damping parameter. The Gilbert 
damping was introduced originally on a phenomenological basis. However, in a conducting ferromagnet system, the 
damping parameter can be derived microscopically [19]. In the case of a ferromagnetic BEG, the damping term arises 
due to collision with non-condensed atoms, which is introduced in a phenomenological manner. 

The analogy between the dissipative hydrodynamic equation and the extended LLG equation implies interesting 
connections between ferromagnetic BEGs and conducting ferromagnets. For instance, the current-driven motion of 
domain walls and spin vortices, which has been investigated in conducting ferromagnets theoretically [17, 20, 21] and 
experimentally [22, 23], can be investigated also in ferromagnetic BEGs by comparison. More interesting phenomena 
such as the anomalous Hall effect, which is the Hall effect due to the magnetization in a conducting ferromagnet [24- 
27], may be investigated in a ferromagnetic BEG from the viewpoint of the interaction between current and spin 
configuration. Since there are no impurities in a ferromagnetic BEG, which is also indicated by the absence of the 
nonadiabatic term in the; hydrodynamic equation, wc; can expcict to investigate pure adiabatic spin-transfer effects in 
this system. These interesting connections motivated us to investigate the domain wall motion and the efi'ect of the 
superfluid current in a ferromagnetic BEG. In this paper, we demonstrate the magnetic domain pattern dynamics, 
which arc mainly simulated by the dissipative hydrodynamic eqiiation. 

In the study of magnetization dynamics, we take into account the magnetic dipole-dipole interaction (MDDI) and 
the quadratic Zeeman effect. The MDDI is known to yield the spatial structure of magnetizations [28-31]. On the 
other hand, the quadratic Zeeman energy determines the easy axis of the magnetization. In this paper, wc consider a 
quasi-two dimensional (2D) system and choose the easy axis normal to the 2D plane. Labyrinthine or striped patterns 
then appear in the magnetic domains with the magnetization parallel and anti-parallel to the normal direction, 
similar to the domain patterns in a ferromagnetic thin film [32]. We numerically investigate the domain formation 
dynamics with and without the superfiuid current, and find that the superfluid current helps the spin transport to 
reach a stationary configuration. We also discuss characteristic lengths of a domain pattern in the stationary state. 
The analytical estimation of domain size shows good agreement with the averaged domain size of the numerically 
simulated domain pattern. 

The paper is organized as follows. In Sec. II, wc derive the dissipative hydrodynamic equation from the GP equation 
with dissipation. In Sec. HI, we explain how the MDDI is implemented in the hydrodynamic equation and which 
type of magnetic domain pattern is expected to appear depending on the balance between the MDDI energy and 
the quadratic Zeeman energy. We demonstrate the dynamics of the magnetic domain patterns simulated by the 
dissipative hydrodynamic equation in Sec. IV. The time evolution of average longitudinal magnetization, kinetic and 
MDDI energies is also shown for different quadratic Zeeman energies. The domain pattern dynamics are compared 
between hydrodynamic equation simulations with and without the superfiuid current and those of the GP equation. 
In Sec. V, we theoretically estimate the characteristic lengths of magnetic domain patterns. The characteristic domain 
size shortly after the emergence of a pattern is estimated from the dynamical instability, and that of the stationary 
pattern is estimated using an ansatz of a stripe domain configuration. Gonclusions and outlook are given in Sec. VI. 



Wc consider a spin-i^ BEG of TV atoms under a uniform magnetic field applied in the z direction confined in a 
spin- independent optical trap i7trap('")- The zero-temperature mean-field energy is given by 



where Sy^im Strap, £p, £-q, £-s, and E^a are the kinetic energy, the trapping potential energy, the linear and quadratic 
Zeeman energies, the short-range interaction energy, and the MDDI energy, respectively. The kinetic and the trapping 
potential energies are given by 



II. DISSIPATIVE HYDRODYNAMIC EQUATION 



^ — ^kin H" ^trap H" ~^ ^ddj 



(3) 




--F 



(4) 




F 



(5) 



--F 
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respectively, where ^'m(r) is the condensate wave function for the atoms in the magnetic sublevel m. The wave 
function is normahzed to satisfy 



= fdr f2 (6) 



N ^ 

-F 

The hnear and quadratic Zeeman energies under the external magnetic field B = Bz are given by 

f F 

£p=p dr Yl (7) 

•' m,n=-F 

r ^ 

e, = q dr '^*ra{r){F^.)mn'^n{rl (8) 

^ m,n=-F 

respectively. Here, Fx,y,z arc the spin-_F matrices. The linear Zeeman energy per atom is given by p = gpfi^B, where 
qf is the hyperfine ^-factor, and /Ub is the Bohr magneton. The quadratic Zeeman energy is induced by a linearly 
polarized microwave field as well as by an external magnetic field: q = gs + 9em, where qb = (fl'FMB-B)^/-Ehf , with Eh[ 
being the hyperfine splitting energy, and ^em = — fv^^"^/ (4^), with f2 being the Rabi frequency and 5 the detuning [33]. 
The short-range interaction energy is 

^« = 2 / '^''5^ E E E -^as{Fm,Fm'\SMs){SMs\Fn',Fn)^r.'{r)^r,{r), (9) 

m,nm',n' S=0,even Ms = — S 

which comes from the short-range part of the two-body interaction given by 

2F ,2 

V,{ry) = S{r-r') Yl ^""sTs, (10) 

S=0,even 

where Vs = X]ms=-S' \^'^s){SMs\ projects a pair of spin-F atoms onto the state with total spin S*, as is the s- 
wave scattering length for the corresponding spin channel S, and {Fm,Fn\SMs) in Eq. (10) is the Clebsch-Gordan 
coefficient. 

The MDDI energy is given by 

SA6.= ^-f f drdr' Y Ur)Q^,{r-r')Mr'), (11) 

where Cdd = Mo(5FMB)^/(47r), with /io being the magnetic permeability of the vacuum, Qi^v{r) is the dipole kernel, 
whose detailed form is given in the next section, and 

F 

U{r)= Y *™(^)(^M)mn*„(r) (12) 

m,n=— F 

is the spin density. The number density is defined by 

F 

ntot(r)= Y (13) 

m=-F 

The hydrodynamic equation without dissipation has been derived from the GP equation [6-8]. Here, we consider 
dissipation, which can be phenomenologically introduced to the GP equation by replacing id/dt with {i — V)d/dt [34]. 
The origin of the dissipation can be interpreted as the relaxation process of the thermal particles into the conden- 
sate [34, 35]. The value of F is often taken to be 0.03, and, in fact, experimental results have been well explained by 
the dissipative equation with F = 0.03 [34, 35]. The dissipative GP equation is given by 

(.-r)4*.(.,*) = ^(^-^'^W) 



E 

n=-F 



*n(r,i), (14) 
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where we use the time-dependent chemical potential so that the total number of atoms is conserved. The 
spin-dependent part i?„i„ is given by 



Hmn{r,t) = Utri,p{r)S,nn + p{Fz)rnn + ^(^f )r 



F 2F S 

E 

■m',n' = -F S=0,ovon Ms = -S 



E E 

' S- 
+ Cdd 



as {Fm, Fm'\SMs) {SMs\Fn', Fn) (r , (r , t) 



(15) 



where the non-local dipole field 6(r, t) is defined by 



Q^Ar-r')Mr',t). 



(16) 



Below we omit the summation symbol; GrcH-k indices that appear twice are to be summed over x, y, and z, and 
Roman indices are to be summed over —F,...,F. 

In the following, we consider a ferromagnetic BEC. We assume that the BEC is fully magnetized, |/| = Fritot, and 
only the direction of the spin density can vary in space. This assumption is valid when the ferromagnetic interaction 
energy is sufficiently large in comparison with the other spinor interaction energies, MDDI energy, quadratic Zeeman 
energy, and the kinetic energy arising from the spacial variation of the direction of /. The linear Zeeman effect is not 
necessarily weaker than the ferromagnetic interaction, since it merely induces the Larmor precession. For example, 
the short-range interaction (10) for a spin-1 BEC can be written as (see Appendix B2) 



{mn\Vs{r, r')\m'n' 



where cq — 47r?i^(2a2 -I- ao)/{3M) and ci 



= 5{r 



r') [coSr, 



(17) 



02 



ao)/{3M). The ground state is ferromagnetic for ci < 0. 



The above assumption is valid when \q 
larger than the spin healing length ^sp 



|ci|ntot; Cdd |ci|, and the length scale of the spatial spin structure is 



h/ y^2M\ci\ntot- Moreover, in the incompressible limit, namely when the 
spin-independent interaction (cgritot for the case of a spin-1 BEC) is much stronger than the ferromagnetic interaction 
and MDDI, the number density Utot is determined, regardless of the spin structure, and assumed to be stationary. 
This is the case for the spin-1 ^^Rb BEC. 

We introduce a normalized spinor with '^/jn{r,t) — \/ntot{r, i)Crn('") t)- When the atomic spin is polarized in 



the z direction, the order parameter is given by Q 
the gauge transformation and Euler rotation as 



(0) 



SmF- The general order parameter is obtained by performing 



Cr. 



J(4>-Fy)^-tF,a^-tFyl3^(0) 



(18) 



where a, (3, and 7 are Euler angles and is the overall phase. Due to the spin-gauge symmetry of the ferromagnetic 
BEC (i.e., the equivalence between the phase change ^ and spin rotation 7), distinct configurations of Cm ai'e charac- 
terized with a set of parameters a, /?, and cj)' = <j) — Fj. The unit vector of the spin density, / = f / {Futot), for the 
order parameter (18) is denoted by a and /3: 



(0) 



''sin /3 cos a^ 
sin ^ sin a 
cosP 



where TZ is an S0(3) rotation matrix given by 



'cos a —sin Of 0' 
71= I sin a cos a 
1, 



cos/3 sin/3^ 

1 
— sin/3 cos/3. 



(19) 



(20) 
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The time evolutions of the total number density and the normalized spin density are given by 



dnt 



dt 

df 1 



d_ 
dt 



\1/ 



dt Frit 



d_ 

dt 

d_ 

dt 



dt 



V]/* 



m I ^ run ^ n 



dntot f 
dt ntot ' 



respectively. The velocities of the superfluid current and spin superfluid current are defined by 

h 



:[C(Va) - (VC)Cn 



V 



'^P™ 2Mi 



2Mi' 

A[v</,'-F(Va)cos/3], 

h 



(^^MU[C(vCn)-(vc)C. 



respectively. Substituting Eq.(14) into Eq. (21), we obtain 

dnto% 1 



where 



dt 



Mlocal 



1 + r 



rV • (ntot^mass) 



r 2 



1 + 



— ir 

2 



ntot 2M ^ntot 

From the above equations, we can also derive the time derivative of v- 



+ — <ass + ^(V/)'. 



dt 



1 + r2 2ntc 



-V • (ntot^^mass) 



+ hF{S7f)-lfx^ 



The detailed derivation is given in Appendix A. 

Next, we consider the incompressible limit and assume dntot /dt = 0. Then, Eq. (25) leads to 



V • (ntott'mass) = -^rntot [/^local " l^{t)] ■ 



This equation simplifies Eq. (27): 



d_ _ 

dt 2ntotr 



^ V [V • (ntot f mass)] + hF(yf) • ( / X ^ ) • 



Substituting Eq.(14) into Eq. (22) and using Eq. (28), we obtain the equation of motion for spin as follows: 



df _ 1 

dt 1 + r2 



-f X Beff - {V„ 



■V)/ 



r 



i + r2 



-/ X Beff - {Vn 



V)/ 



:(a • V)/ - Cddb -pz- q{2F - l)fj, 



(21) 
(22) 



(23) 



(24) 



(25) 



(26) 



(27) 



(28) 
(29) 

(30) 
(31) 



2M 2M' 

where a = (Vntot)/ntot- The detailed derivation is given in Appendix B. The effective field Bgff is also derived from 
the reduced Hamiltonian 'Hmag that contains only spin-dependent terms: FB^^ = —SHma.g/Sf, where 



"^mag 



1 

ntot 



dr 



(V/)" , Cdd^ , ^ , g (2F-l) ,„ 2 



(32) 



and b includes /. 

Here we note that Eq. (30) is equivalent to Eq. (2) and corresponds to the extended LLG equation, which describes 
the magnetization dynamics in a conducting ferromagnet in the presence of spin currents interacting with magnetiza- 
tion, with the adiabatic spin-transfer torque. The superfluid velocity t'mass in Eq. (30) corresponds to the spin current, 
which is associated with the electric current density in the extended LLG equation of a conducting ferromagnet. 
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III. MAGNETIC DIPOLE-DIPOLE INTERACTION 



Below, we consider a BEC confined in a quasi-2D trap whose Thomas-Fermi radius in the z direction is smaller 

than the spin healing length. Wc approximate the wave function in the z direction by a Gaussian with width d: 
*rTi(^-L)-2^) = V'm(''x)/j(-z)) where rj_ = {x,y), and h{z) = exp[—z^/{4:(P)]/{2TTcPy^^. When we consider a quasi-2D 
BEC, ntot, /, and Umass are defined by means of 'tl)m instead of [36]. If one replaces with 'tl)m,, as with rjas, 
Cdd with r/Cdd, and b with 5, the equation is the same as Eq. (14), where rj = J dzh'^{z)/ f dzh?{z) = l/\/47rd^ and 

K{r±) = I dVlO(fD)(r^ - rl) [rm{r±KF.)mnMr'±)] , (33) 

with 

Q^^^\r^ -r'^) = l^JJ dzdz'h\z)h\z')Q^,{r - r'). (34) 
The 2D dipole kernel in the laboratory frame of reference is given by 

g(2D,Iab)(^^^^^....g(2D,ab)^ (35) 
k 

where the subscript _L is omitted for simplicity and 

^^2D,iab) ^ _^ I ^ \+47rG{kd)\ kjy P I, (36) 




with k = {kx, ky), k = \k\, k^.y = k^.y/k, and G{k) = 2ke^ e~* dt = ^/nke'^ erfc(fc). It can be shown that G{k) 
is a monotonically increasing function that satisfies G(0) = and 6(00) = 1. 

When the linear Zeeman energy is much larger than the MDDI energy, we choose the rotating frame of reference in 
spin space by replacing with e~*^'"*/'*^'m, and eliminate the linear Zeeman term from the GP equation. In this 
case, the contribution of the MDDI is time averaged due to the Larmor precession. The 2D dipole kernel, which is 
averaged over the Larmor precession period, under an external magnetic field in the z direction, is given by [36] 

Q(2D,rot) ^ (^^^ _ 3^^^^^^) ^ ^^k■r■ q^^ (37^ 

k 

where 

Qfe = — [-2 + 3G(fcd)]. (38) 

The formation of a stable magnetic domain pattern depends on the quadratic energy as well as the MDDI. Since 
the quadratic Zeeman energy is the monotonic function of q{z ■ f)^, transverse {fz = 0) and longitudinal {fz = 
±1) magnetization is preferable for g > and q < 0, respectively. For <; > 0, the uniform pattern of transverse 
magnetization is stable because of the MDDI. When a strong magnetic field is applied, the dipole kernel is given by 
Eq. (37), which is low for small-fc modes of a transverse magnetization pattern and for large-fc modes of a longitudinal 
one. Thus, a imiform transverse magnetization pattern (i.e., fc = 0) is stable for q > 0. The situation is the same for 
a magnetic pattern under zero field, in which the dipole kernel is given by Eq. (35). For g < 0, the uniform transverse 
magnetization pattern is still stable if the quadratic Zeeman energy is smaller than the MDDI energy. When q < 
and \q\ is large enough to balance with the MDDI, the uniform transverse magnetization becomes unstable and a 
longitudinal magnetization pattern appears. In other words, there is a threshold of q where a non-uniform longitudinal 
magnetization pattern appears. The threshold can be calculated from the linear stability analysis, as will be discussed 
in Sec. V. 

Here, we focus on the negative-q regime, which can be achieved by means of a linearly polarized microwave field even 
in the absence of an external magnetic field [see, below Eq. (8)]. When g < and \q\ is sufficiently large, longitudinal 
magnetization is dominant and magnetic domains with ~ ±1 form patterns because of the MDDI. This situation 
is consistent with that of a uni-axial ferromagnet [32] . 
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FIG. 1; (Color online) (a) Time dependence of the averaged longitudinal magnetization, kinetic and MDDI energies in the 
presence (solid curves) and absence (dashed curves) of Vmass- The dot-dashed curve in the graph of E-^i^/h (middle panel) 
corresponds to the contribution from Vmass- (b) Snapshots of longitudinal magnetization in which white and black correspond 
to positive and negative values of fz, respectively. The size of each snapshot is 256 x 256 /^m. Here, we take q/h = —50 Hz, 

and ntot = ^^uli where ri'^J^ = y/2TTd?n'^^ with — 2.3 x 10^** cm~'^ and d = 1.0 ^m. The damping rate is given by a typical 
value F — 0.03. The other parameters are given by the typical values for a spin-1 *^Rb atom: M = 1.44 x 10"^^ Kg, F = 1, 
and qf = —1/2. 



IV. DOMAIN FORMATION DYNAMICS 



In this section, we focus on the domain formation dynamics with and without Vmass under a strong magnetic field 
to illustrate how fmass affects domain pattern formation. Later, the domain patterns under zero magnetic field are 
also shown as complementary results. The validity of the dissipative hydrodynamic equation is also examined by 
comparing hydrodynamic and GP equation simulations. 

For the magnetic domain pattern simulations, we solve the coupled equations (29) and (30), where Vmass and Bcs 
are defined by Eqs. (23), and (31), respectively. For simplicity, we take a = and use periodic boundary conditions. 
In the case of a strong magnetic field, we employ Eq. (37) as a dipole kernel. To see the role of t'mass, we demonstrate 
the calculation without Vmass, in which Vmass is always taken to be zero, as well as the full calculation using all those 
equations in the presence of Wmass- The initial condition is fx \^ fy ^ 0, and ~ with small noises, and t^mass = 0. 

Figure 1(a) shows that the longitudinal magnetization grows rapidly to form magnetic domains. The kinetic and 
MDDI energies also grow rapidly at first. After the rapid increase, the kinetic energy decays as magnetic domains grow 
and the domain wall density decreases. The averaged longitudinal magnetization, kinetic energy and its contribution 
from Vmass , and MDDI energy (per unit area per atom) are defined by 



dd 



1 
l2 



d'r |/,(r)|. 



1 

l2 



1 
Z2 



^2 r 



AM 



1^ Y^mass(r) 



(V/.(r))2 + (V/,(r))2 + (V/.(r))^ 



(fr ^?76(r) • /(r). 



(39) 
(40) 
(41) 
(42) 



respectively. Here, L is the system size. 

The longitudinal magnetization grows rapidly for a short time, and then the averaged longitudinal magnetization 

saturates at around \fz\ — 1 [see the top panel of Fig. 1(a)]. Figure 1(b) shows that magnetic domains emerge in a 
short time and, after that, spread slowly to reach a stationary configuration. After the rapid growth of fz, magnetic 
domain patterns grow faster in the presence than in the absence of Vmass- The difference is apparent in the series 
of snapshots and -Ekin and Edd, although the contribution of t'mass to -Ekin is small compared with the total kinetic 
energy [see the middle panel of Fig. 1(a)]. In the presence of Vmass, magnetic domains grow efficiently because of spin 
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FIG. 2: (Color online) Comparison of domain formation dynamics between the hydrodynamic description and CP numerics, 
(a) Time dependence of the averaged longitudinal magnetization, kinetic and MDDI energies simulated by the hydrodynamic 
equation (left column) and the GP equation (right column), (b) Snapshots at f = 5 s simulated by the hydrodynamic and GP 
equations. Solid, dashed, and dot-dashed curves in (a) are for q/h = —50, —40, and —30, respectively. 



transfer accompanied by the transfer of atoms by means of i^mass ■ 

To examine the vahdity of the hydrodynamic description, we compare the simulations of the dissipative hydrody- 
namic and the GP equations. In both cases, the averaged \fz\i E^im E^di and snapshots for several values of q are 
shown in Fig. 2. For the GP equation, the averaged longitudinal magnetization, kinetic and MDDI energies (per unit 
area per atom) are defined by 













i.l(_. 






Edd = 





\h{r)\ 
ntot{r) 

2M J 

'\ ^Vb{r) ■ f{r), (45) 



(43) 

dV^<„(r)VV™(r), (44) 



respectively. Here / and ntot are defined in Eqs. (12) and (13), respectively. These values are equal to Eqs. (39), (40), 
and (42) if ritot is uniform and |/(t')| = fitot- Below the threshold where a non-uniform pattern begins to appear, 
the most stable pattern is a uniform pattern. The threshold is q/h ~ —30 for the hydrodynamic description and 
q/h ~ —40 for the GP equation. Near the threshold, the averaged longitudinal magnetization is small compared to 
those in other cases. Except for in the vicinity of the threshold, the time dependence of domain pattern formation looks 
similar for the two simulations. However, the domain size differs significantly due to differences in the domain wall 
structure. In the GP simulation, the amplitude of magnetization is suppressed in domain walls, while the suppression 
of magnetization does not occur in the hydrodynamic description. An analysis of domain sizes for both cases will be 
given in the next section. 

The magnetic domain patterns under zero field are simulated by the dissipative hydrodynamic equations with 
Eq. (36). In Fig. 3, domain pattern formation from the initial condition of ~ 1 is shown in the presence and in 
the absence of Vmass- The magnetic domain patterns under zero field look similar to those under a strong magnetic 
field. The effect of t>mass is also similar: Vmass moves the domain walls faster. However, they are strongly affected by 
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time (sec.) 0.1 0.5 1.0 2.0 



no V 




FIG. 3: Domain pattern formation in the absence of a magnetic field. The direction of the stripes depends on the initial 
configuration of spins. The parameters are the same as those given in Fig. 1 



the initial condition. This is because the spin and orbit degrees of freedom couple in the MDDI under zero field [see 
Eq. (36)]. If the initial condition is given by /j, ~ 1, one can see remarkably similar domain patterns as shown in the 
snapshots in Fig. 3, which are rotated by 90°. 



V. CHARACTERISTIC LENGTHS OF DOMAIN PATTERNS 



We have seen in the previous section that the magnetic domain pattern has an initially short characteristic length 
that later increases in size. At the beginning of domain pattern formation, domain size is estimated from the linear 
instability analysis. The domain size of the stationary pattern can be estimated as that of a stripe domain pattern. 
For the parameters given in our numerical simulations, the domain size of the stable stripe pattern is longer than the 
domain size estimated from the linear instability. 



A. Dynamical instability 



The dynamical instability (linear instability) under a strong magnetic field has been discussed previously, both for 
the hydrodynamic equation [7] and for the GP equation [36] . For the hydrodynamic equation, the growth rate of the 
unstable mode is calculated from the eigenvalues of the linearized equation of small deviations from the uniform initial 
condition. For the GP equation, a similar equation is derived by means of Bogoliubov analysis. When the initial 
magnetic pattern is uniform and fx = 1, the respective growth rates of the unstable mode for the hydrodynamic and 
GP equations are given by 




q - 47rcdd"'tot [1 - G{kd)] 



2M 



2'KCddntotG{kd) 



h^k^ 



q - 47rCddntot(l - a)[l - G{kd)] 



h^k^ 



2M 



+ 27rCddntot(l + q)G{kd) 



where fitot = ntotf] and 



q = 



2ritot(|ci| + 47rcdd/3) ' 



Here, we have considered Wmass = and neglected the dissipation, which is very small (i.e., T <C 1). 
The domain size at the emergence of a non-uniform pattern is estimated as 

where ko is given by the momentum at which A^(A;o) or A'^(feo) has its maximum value. 



(46) 
(47) 

(48) 

(49) 



B. Domain size estimated from the hydrodynamic equation 



Here, we estimate characteristic lengths of the domain pattern in the stationary state. For q < with large \q\, 
the ideal stable pattern is the stripe pattern of longitudinal magnetization. We assume that the stable pattern is 
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described by 

= cn(a;/K^, k?), fy = 0, f^= sn(a;/«;^, k?), (50) 

where sn(a;/fi;^, K^) and cn(.i;/K^, k?) are the Jacobi ehiptic functions with < < 1. These functions contain the 
characteristic lengths of the stripe domain pattern [37]: the domain wall width ^ and the periodicity of the pattern 
24 with 

4 = 2k^K{k'^), (51) 

where K{k'^) is the complete elliptic integral of the first kind. The kinetic and quadratic Zeeman energies are calculated 
as 

^'''"=2Mi4^' ^^^^ 



(53) 



where E{k^) is the complete elliptic integral of the second kind. The domain wall width is estimated to be the length 
at which the summation of the two energies has a minimum value. Solving 9(£'kin + Eq)/d^ = 0, we obtain 

For estimation of 4, we need to take the MDDI energy into account. For simplicity, we assume 4 ^ and take 
fx = fy = and 



f _ / -1 foi' (2" - 1)4 < a; < 2n4 
■'^ ~ \ 1 for 2n4 < a; < (2n + 1)4 



OO 



1 _ /I ^n 



where n is an integer. From Eqs. (33), (37), and (42), we can calculate the MDDI energy as 

-(-ir " 



-E'dd = — ^Cdd^tot [-2 + 3G{mTd/£s 

n=l 



mr 



(56) 



Instead of using the original definition of G{k), we use an approximate function G{k) ~ 1 — exp(— ^/ttA;). Then, 
Eq. (56) is rewritten as 



47r , ri 6 
■C'dd = --;rCdd"-tol s t: H — 7 



Li2(-e-'^^'^/^0 -Li2(e-'^^''/^'')l I, (57) 



where Lis(z) = ^2^=0 ■2'^/^* is the polylogarithm. We estimate 4 to be the length at which the total energy has its 
minimum value in the limit of k — )• 1 [thus E{k'^) 1]. Solving 9(i?kin + Eg + Edd)/dis = and taking k = 1, we 
obtain 

4 = ^ . / (58) 



arctanh 



exp ( ^IlMi^] 

\ ScddfitotV^d J 



In Fig. 4(a), we plot the q dependence of 4,C: and £s for ntot = Qn['^^. The estimated domain size is in good 
agreement with the average domain size of the simulations; e.g., for q/h = —40, the estimated domain size £s ~ 24 
/Km from Eq. (58) and the average domain size 4 ^ 20 nm from simulations shown in Fig. 2(b). The stable domain 
size 4 is larger than the initial domain size £i. This property is consistent with the simulations in which initially 
small magnetic domains appear within a short time before spreading to form a stable (or metastable) pattern. The q 
dependence of 4 is also consistent with the simulations in which the domain size increases as jgj increases. 

The domain size of a stable pattern depends also on the number density. In Fig. 4(b), 4 is plotted for n^ot = '^f^ioL 
6n^Qj, lOn^Qj and 20?!^^^. The smaller the value of ntot, the larger the £s. For the typical number density ritot = ''^tot 
in experiments [3-5], the estimated domain size is too large to observe in a conventional experimental system. 
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20 40 60 



-q/h 



FIG. 4: (Color online) Dependence of theoretical characteristic lengths on ~q/h for the hydrodynamic description, (a) Theo- 
retically estimated characteristic lengths of a stripe pattern (solid curve) [Eq. (58)], a pattern at its emergence (dashed curve) 
[Eq. (49) with Eq. (46)], and the domain wall width (dot-dashed curve) [Eq. (54)] for ntot = Bn^Qt. (b) Theoretical characteristic 

length of a stripe pattern [Eq. (58)] for ntot = 3n{o| (dashed line), Qri^J^ (solid line), lOnj^t (dot-dashed line), and 2QnfJ^ (dotted 
line). 




FIG. 5: (Color online) Dependence of theoretical characteristic lengths on ~q/h for the GP equation, (a) Theoretically 
estimated characteristic lengths of the stripe pattern with fully-magnetized domain walls (solid curve) [Eq. (58)], that of the 
stripe pattern without transverse magnetization (dot-dashed line) [Eq. (65)], and that of a pattern at its emergence (dashed 
curve) [Eq. (49) with Eq. (47)] for ntot = 6n'°'. (b) Theoretically estimated domain width of the stripe pattern with fully- 
magnetized domain walls (solid curve) [Eq. (54)] and that of the stripe pattern without transverse magnetization (dot-dashed 
line) [Eq. (64)] for ntot = 6n^2- 



C. Domain size estimated from the GP equation 

In the GP equation, spins are not always fully magnetized. Even when ~ ±1 over most of a magnetic domain 
pattern, the magnetization in the domain walls may vanish. The domain wall structure depends on the ferromagnetic 
interaction (cintot), which competes with the quadratic Zeeman and kinetic energies. The ferromagnetic interaction 
lowers energy for high spin density. In other words, the transverse magnetization exists in the domain wall between 
magnetic domains with Jz — ±1, when the ferromagnetic interaction is strong. On the other hand, the quadratic 
Zeeman energy prefers a state with sublevels m = ±1. Here, let us consider two states with = 0; for instance, 
(a) (Ci,Co,C-i)^ = (1/2,1/72,1/2)'^ and (b) (1/^2, 0, l/\/2)'^. State (a) corresponds to ^ I and fy = = 
(transverse magnetization), and state (b) to fx = fy = fz = (zero magnetization). Comparing them, one can find 
that the quadratic Zeeman energy of state (b) is smaller than that of state (a). In other words, domain walls with 
zero magnetization appear when the quadratic Zeeman energy is dominant. Therefore, the domain wall structure 
changes at around |ci|ntot — \q\- 
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Now we estimate characteristic lengths of magnetic domains for |ci|ntot < I?]- Since zero magnetization is preferred 
in domain walls, we assume the stripe pattern in this case is described by 



/ritot cos ■ 



7r/2 



ipo = 0, V-i 



/ntot sm 



2 , -r^ -r . V— 2 ' 

with 9 = am(a;/K^, K^), where a,m.{x /k^,k^) is the Jacobi amplitude. Equation (59) corresponds to 

fx = fy = 0, fz = sn(x/<, K^). 
Then, the kinetic, quadratic Zeeman, and ferromagnetic interaction energies are calculated as 

1 E{k^) 



kin 



Eg = q, 



cl 



cintot 1 
2 L2 



/ 



dr'' f{r) 



cintot 



K 



1 2g ^(k^) 
"2 



Note that is independent of ^ and 4- Solving 9(£'kin + Eci)/d^ = in the «; ^ 1 limit, we obtain 



A/-4Mcintot " 



(59) 

(60) 

(61) 
(62) 

(63) 
(64) 



In Fig. 5(b), Eq. (64) is plotted as the dot-dashed line, while the solid curve is the plot of Eq. (54). Equation (64) is 
vahd for |g| > |ci|fitot — 35/i. In this region, the values of ^ of both equations almost coincide, as seen in Fig. 5(b). 

To estimate the domain size of the stripe pattern, we again assume the MDDI energy is given by Eq. (57). Solving 
9(Skin + Ed + Edd)/d£s = with k = 1, we obtain 



4 = 



1 



arctanh 



exn ( 

^\ 16cddratot\/7r<i/ 



(65) 



In Fig. 5(a). Eq. (65) is plotted as the dot-dashed line and compared with the solid curve given by Eq. (58). The 
dashed c;urve expresses the domain size £i at the beginning of domain pattern formation and has a nonzero value 
above the threshold. Since l^] > |ci|ntot in the region above the threshold, the domain size 4 is estimated by Eq. (65) 
instead of Eq. (58) for the stable pattern simulated by the GP equation. Also in this case, the estimated domain size 
is in good agreement with the average domain size obtained from simulations. 



VI. CONCLUSIONS AND OUTLOOK 



We have derived the dissipative hydrodynamic equation of a ferromagnetic Bose-Einstein condensate (BEG). This 
equation has the same form as the extended Landau-Lifshitz-Gilbert (LLG) equation, which was originally developed 
to explain the spin dynamics in a conducting ferromagnet interacting with spin-polarized currents, including an 
adiabatic spin-transfer torque term. The dissipative hydrodynamic equation enables us to investigate how the domain 
formation dynamics are affected by the superfluid velocity Umass, which is inseparable in the Gross-Pitaevskii (GP) 
equation. 

We have demonstrated domain pattern formation simulated by the dissipative hydrodynamic equation with and 
without t^mass- Although no remarkable difference appears at the beginning of domain pattern formation, Vmass has 
an effect on later domain formation dynamics: pattern formation is faster in the presence than in the absence of f mass- 
We have also shown simulations of the GP cqiiation and compared them with those of the hydrodynamic equation. 
The dependence on q, which characterizes the quadratic Zeeman energy, of domain pattern formation is different 
between hydrodynamic and GP simulations because the threshold for a nonuniform magnetic pattern differs between 
them. Nevertheless, for large \q\, magnetic domain patterns eventually come to look similar in both simulations, 
although they differ in size. The difference is caused by the fact that the assumption of full magnetization for the 
hydrodynamic description is not satisfied for large \q\ in the GP simulations. 
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To explain the difference in domain size between the hydrodynamic and GP equation simulations, we have estimated 
the characteristic lengths of the domain patterns. The domain size at the beginning of pattern formation is estimated 
by means of the linear stability analysis. The difference in domain size for a short time is based on the dynamical 
instability. The domain size of the domain pattern in the stationary configuration is estimated using the ansatz of 
the stripe pattern of longitudinal magnetization domains. The analytical estimations are in good agreement with the 
numerical simulations. 

In conclusion, the dissipative hydrodynamic equation provides a simple approach to discuss magnetization dynamics 
in a ferromagnetic BEG. The hydrodynamic equation simulation qualitatively well reproduces the domain formation 
dynamics that are simulated by the GP equation. However, quantitative discrepancies arise when the assumption 
of the hydrodynamic description fails: for instance, the ferromagnetic interaction energy becomes comparable with 
the quadratic Zeeman energy. The analogy between the dissipative hydrodynamic equation and the extended LLG 
equation can provide suggestions on new experiments of a ferromagnetic BEG to investigate interesting phenomena 
that are observed in conducting ferromagnets. 

One of such interesting and possible phenomena is the anomalous Hall effect (AHE), which is the Hall effect due 
to the magnetization and observed in conducting ferromagnets. Here, we focus on the AHE caused by a skyrmion 
configuration or spin chirality [24-27], although there are several mechanisms to cause the phenomenon. The Berry 
phase, which is generated by a skyrmion configuration, induces an effective magnetic field or gauge flux. When an 
electric field is applied, electrons move to the perpendicular direction to both the electric and the effective magnetic 
field. This is the mechanism of the AHE due to spin chirality. Inversely, the electric current can move the skyrmion. 
The AHE is expected to be observed very clearly in the adiabatic limit, where a ferromagnetic BEG is supposed 
to realize. Unfortunately, there are difficulties for observation of the AHE in a ferromagnetic BEG: for example, 
it is difficult to create the external field that corresponds to the electric field of a conducting ferromagnet system. 
However, the essential aspect of the AHE (i.e., the interaction between current and spin chirality) can be investigated 
in a ferromagnetic BEG. Experimentally, one can create a current by sudden change of a trapping potential or the 
fictitious field that is created through the vector potential induced by a laser field [38]. The investigation about the 
interaction between current and spin configuration in a ferromagnetic BEG will give an insight on pure adiabatic 
spin-transfer effects. 
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Appendix A: Time evolution of v, 



mass 



Substituting Eq. (18) into Eq. (14) and applying from the left, we have 




(Al) 



where we have used '^m = e'"^ Umn'^n^ ■ The first term on the right hand side is calculated as 




(A2) 



where we have used Eq. (23) and the following equations: 

f,iFy/3p^^-iF,l3 ^ ^.Qg /3 - sin;3, 
F 



C> {[F.(Va) sin^ - Fy{VP)f] = ^ [(Va)^ sin^ /? + (V/?)^ 



The imaginary part of Eq. (Al) reads 

m 



2rn, 



tot 



= ? (V/? 



— - i^COS/3— - V • (ntot'Wmass)- 



Prom the real part of Eq. (Al), we have 

?ir dntot 
2ntot dt 



+ n 



dt 
'dt 



dt J 



F COS 13— ] + [/Xiocal - mW] = 0, 



where 



MlocaKT",*) 



h^F 



ntot 



Combining Eqs. (A6) and (A7) and eliminating {d<p'/dt — F cos (3da/dt), we obtain Eq. (25) as follows: 

dntot 1 



rV • (ntot'Wmass) 



On the other hand, substituting Eq. (A6) into Eq. (A9), we have 



r 2 

T"tot[Mlocal - l^{t)]- 



h 



,da 



h 



— - Fcos/3-j = -^Kcai - m + TTr^^V • {ntotv. 
The gradient of the left hand side of Eq. (AlO) is 



da 



dt 



dt 



^V ^-i^cos^— =h{-{\/<j>' -F{Va)cos(3]-F 



d 



dt' 



d 

M — Vma^ss - hF 

dt 



^,da _ d cos 13 
(Vcos/5)--(Va)^ 



dB da 
(Va)sin/3^-(V^)sin^^ 



Here, we introduce vectors rfi and n, which arc orthogonal to /, 

^ sin /? cos a\ / cos /3 cos 

sin /3 sin a , 771= cos /3 sin a 
cos P / \ — sin ^ 



sm a ' 
n = I cos a 



One can easily see 



Employing these vectors, we have 



m = n X f , n ^ f X rh 
df = {dj3)rh + {da) sin/3n. 



(V/) • / X 



df 
dt 



[(V^)m + (Va) sin I3h] 



5/3 . da 
-m+-sm/3n 



' dj3 da 
[(V/3)Tn + (Va) sin /3n] • ( -^ri — ^ sin /3m 



,9/3 



dt 
da 



dt 



= (Va)sin/3^-(V/3)sin^^^ 
From Eqs. (AlO), (All) and (A15), we obtain Eq. (27) as follows: 

Miocai - i^{t) r h 



d_ 

di' 



M — Wmass = -V 



i + r2 



1 + r2 2nto 



V • (ntot^^mass) 



+ m'^f) {fx 



di 
dt 



Appendix B: Time evolution of / 



The time evolution of the normalized spin density is calculated as 



1 



dt Fntot 
1 



dt 



dntot fa 



dt ritot 



Fntot HI + t 2M 

+2 [Im {^*^{F^)mlHl„^n) - TRe {^*^{F^)mlHin^n)] + 2Tf,{t)} , 



where we have assumed dntot /dt = 0- 



1. Kinetic energy terms 

First, we calculate the spacial derivative terms. Making use of Eqs. (A4) and (A12), we have 

K(^/.)mnV'*„ = Ci°^*(F,m^ + Fyfi^ + fJ^) { V^V^V^ 

+ mtotV • [V0' - {Va)F, cos/3 + (Va)F, sin/3 - iV(3)Fy] 
+ iVntot • [V(/.' - (Va)F^ cos /3 + (Va)F^ sin p - {Vf3)Fy] 

-ntot [V</)' - (Va)F, cos/3 + {Va)F, sin/3 - {Vl3)Fyf} Ci°^. 
Using Eqs. (23), (A5), (A13), (A14) and the relation 

(Va)sin/3m - (V/3)n = -/ x V/, 
(Va) sin/3n + (V/3)m = -/ x (/ x V/), 



we obtain 



Fn, 



tot 



V^-lot 



M 



+ i- 
i 

~ 2 



.M 



f x (/ X VV") 

V • (ntot^^mass 



1 



/X (/x^-V/ 

ntot 



"■tot 



/xXl^.V/ 

ntot 



From this equation and Eq. (28), wc can rewrite the space derivative terms of Eq. (Bl) as 



Re 



Im 



Fnt 



2M f^KJ^rnn^ _ , , 

[ ntot WocalJ/M 



M 



-Fntot 

where a = Vntot/ntot- 



(l^mass • V)/ - y X (VV^ + (a • V)/)] I , 

(fmass • V)/ - y X (V^/ + (a • V)/ 
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2. Short-range interaction terms 

Next, we calculate the contributions of the short-range interaction ^^{F^)miHi^^ni where 

= E P^'\'^s\Fn', Fn). (B9) 

s=o 

Now we rewrite the matrix elements of the short-range interaction in terms of the spin matrices. The projection 

operator Vs satisfies the completeness relation Ylis'^s = !> that is, 

2F 

E {Fm,Fm'\Vs\Fn',Fn). (BIO) 

<S=0,even 

On the other hand, from the identity equation 

fe 



{Fi ■ F^f = 
we obtain 

{Fi/^Fj^^ • • • -^i/;, )mn(-^i/i -^1^2 * * * Fit^^m'Yi' 
5=0, even 

Using Eqs. (BIO) and (B12), can be generally expressed as 



(Bll) 



{Fm,Fm'\Ps\Fn',Fn). (B12) 



fe=l 

where Ag and are given by the hnear combinations of as- For the case of F = 1, for example, we obtain 
Ao = Co = ■nh^{2a2 + ao)/(3M) and Ai = ci = 47r/i^(ao - a2)/(3M), and C^T' can be written in the following form: 

= Co5mn^m'n' + Cl{F)mn ' {F)m'n' ■ (B14) 

Substituting Eq. (B13) to Eq. (B8), we obtain 

^^^„ = c^y*;;^'*n' (bi5) 

F 

= Aontot^mn + '^^kntotMvi„2...v^{F„^F„^ ■ ■ ■ F^^)mn, (B16) 

fe=l 

where 

and we have used = ^ntotCm and Cm Cm = 1- 

When wc consider the ferromagnetic state (i.e., when the order parameter is given by = y/ntote^'^' UmnQn^ with 
Cm'' = <5mF), we have 

= 7^.,.^7^.2.^ • • • t^^.^'MZ,,,,^'^, (bis) 
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where = [Fij^Fi,^ ■ ■ ■ F^^)ff, and TZ is defined in Eq.(20). Then, ^^H^^^n is shown to be independent of 

the local spin direction: 



,(0) 



(B19) 



fe=i 



In a similar manner, we obtain 



fe=i 



Note here that F^ or Fy has to appear an even number of times in the product of F^^F^^ ■ ■ ■ F^^. so that M.'^}v2---i'k — 

(B21) 



[F^^F^^ ■ ■ ■ Fi,^)ff is nonzero. Hence, Mv'iv2---vkJ^^I!i'viU2---Vk becomes nonzero only when jj = z: 



■'^'■VlV2---Vk-'^'-H'viV2---Vk ~ "ll'zJ^ J^'-VlV2---Vk-'^'-VlV2---Vk- 



AO) 



AO) 



Then, Eq. (B20) is written as 



Thus, we obtain 



(0) 

UlV2---'^k 



fe=l 



lm['^>*^{F^W)^r^■^n] = 0. 



(B22) 
(B23) 



(B24) 
(B25) 



3. Other terms 

Finally, we calculate the remaining terms. Using Eqs. (15) and (B25), we easily obtain 



Im 



^;^{F^)mlHln'^n 



Fnt 



1 



2iFntot 
1 



{p[F^, F,] + q[F^, F^] + Cdd6.[F^, F.]}^^ M/„ 



2iFntot 

[/ X (pz + q{2F - l)f,z + Cddb) 



Employing Eqs. (15), (B24), and (B25), we obtain 



Re 



Fntnt 



""^ " /^ + f/trap/^ + f S^, + {2F - 1)/,/^ 
'^tot ^ 



9 
2 

Cdd 

2 



(2f2 - 3F + l)/f + + (2F - l)/,5^, 
6^ + (2i^-l)(6./)/Ml , 



where we have used 



^*m{F^F^ + F.F^)mn^n = Fn,,,5^, + F{2F - l)ntot/^/., 

^*miF^F^ + F^F^)mn^„ = Fritot \{2F - l)fj^, + {2F^ -3F + 1)PJ^ + Ff^ . 



(B26) 



(B27) 



(B28) 

(B29) 
(B30) 
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Equation (B28) is rewritten as 



where 



Re 



Fnt. 



\b* H ^ ~ 1 



"tot 



pz + q{2F - l)f,z + Cddb 
\ {/ • [pz + qC^F - l)f,z + Cddb] }^ 

U-l {/ X [/ X (pz + q{2F - l)hz + Cddb)] } 



"tot 



*;,if„„*„ = ^l,Hl^J>n + ntot {c/trap + pF h + |^ [l + 2(i^ - 



+ Cddi=^& • /} . 



Substituting Eqs. (B6), (B7), (B27), and (B31) into Eq. (Bl), we obtain Eqs. (30) and (31) as follows: 



di 

dt 



1 



i + r2 



-/ X Beff - {V„ 



■V)/ 



i + r2 



-/ X Beff - (W„ 



V)/ 



(B31) 
(B32) 

(B33) 

(B34) 
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